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Lecture 1 


Hyperbolic Conservation Laws 


Much of the material for the first two sets of lecture notes has been taken from a graduate 
course on the numerical solution of conservation laws that I taught at the University of 
Washington in 1988 at at ETH-Ztirich in 1989. The complete text of these notes will be 
published soon by Birkhauser-Verlag in a series of “Lectures at ETH”[50]. I am indebted 
to the people at Birkhauser for allowing me to use some of this material for the present 
lecture notes. I have tried to extract the most important features for this shorter set of 
lectures. For those who are interested, more details can often be found in [50]. 

1.1 Conservation laws 

In this first lecture I will discuss the mathematical structure of hyperbolic systems of con- 
servation laws. A basic understanding of this theory is essential for understanding modern 
numerical methods for the solution of these equations. Several excellent references on the 
mathematical theory of conservation laws exist, including Courant and Friedrichs[22], 
Lax[42], Majda[57], Smoller[66], and Whitham[76]. 

Conservation laws are time-dependent systems of partial differential equations (usually 
nonlinear) with a particularly simple structure. In one space dimension the equations take 
the form 

| U (x,t)+-/(«(x,0) = °- (1-1) 

Here u(x, t) £ lR m is an m-dimensional vector of conserved quantities, or state variables, 
such as mass, momentum, and energy in a fluid dynamics problem. More properly, uj is 
the density function for the jth state variable, with the interpretation that /** Uj(x,t)dx 
is the total quantity of this state variable in the interval [xj,X 2 ] at time t. 

The fact that these state variables are conserved means that Uj(x,t)dx should 
be constant with respect to t. The functions Uj themselves, representing the spatial 
distribution of the state variables at time t, will generally change as time evolves. The 
main assumption underlying (1.1) is that knowing the value of u(x,t) at a given point 
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and time allows us to determine the rate of flow, or flux, of each state variable at (x, t). 
The flux of the jth component is given by some function fj(u(x,t)). The vector-valued 
function /(«) with jth component /,(«) is called the flux function for the system of 
conservation laws. 

1.1.1 Derivation 

To see how conservation laws arise from physical principles, we will begin by deriving the 
equation for conservation of mass in a one-dimensional gas dynamics problem, for example 
flow in a tube where properties of the gas such as density and velocity are assumed to be 
constant across each cross section of the tube. Let x represent the distance along the tube 
and let p(x,t) be the density of the gas at point x and time t. This density is defined in 
such a way that the total mass of gas in any given section from x\ to * 2 , say, is given by 
the integral of the density: 

mass in [xi,X 2 ] at time t = I p(x,t)dx. (1.2) 

Jxi 

If we assume that the walls of the tube are impermeable and that mass is neither created 
nor destroyed, then the mass in this one section can change only because of gas flowing 
across the endpoints xi or x 2 . 

Now let v(x,t) be the velocity of the gas at the point x at time t. Then the rate of 
flow, or flux of gas past this point is given by 

mass flux at (x,t) = p(x,t)v(x,t). (1.3) 

By our comments above, the rate of change of mass in [xi,x 2 ] is given by the difference 
in fluxes at Xi and X 2 : 

d f x 2 

— J p(x,t)dx = p(xi,t)v(x u t)~ p(x 2 ,t)v(x 2 ,t). (1.4) 

This is one integral form of the conservation law. Another form is obtained by integrating 
this in time from ti to < 2 , giving an expression for the mass in [xi,X 2 ] at time <2 > *1 in 
terms of the mass at time and the total (integrated) flux at each boundary during this 
time period: 

f p(x,t 2 )dx = f p(x,ti)dx (1.5) 

Jx 1 Jx\ 

f t2 f * 7 

+ / p(x u t)v(xi,t)dt- / p(x 2 ,t)v(x 2 ,t)dt. 

Jt , Jti 

To derive the differential form of the conservation law, we must now assume that p(x, t ) 
and u(x,t) are differentiable functions. Then using 

r * 2 d 

p(x,t 2 )~ p{x,t x) = J — p(x,t)dt 


( 1 . 6 ) 
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yx2 d 

p(x 2 ,t)v(x 2 ,t) - p(x u t)v(x u t) = J* —(p(x,t)v(x,t))dx (1.7) 

in (1.5) gives 

I/ J 3 + *)”(*>*))} <*xdt = 0. ( 1 . 8 ) 

Since this must hold for any section [xi,x 2 ] and over any time interval [h,t 2 ], we conclude 
that in fact the integrand in (1.8) must be identically zero, i.e., 

Pt + ( pv) x = 0 conservation of mass. (1-9) 

This is the desired differential form of the conservation law for the conservation of mass. 
(Subscripts denote partial derivatives.) 

In general, the equation (1.9) must be solved in conjunction with equations for the 
conservation of momentum and energy: 


{pv) t + (pv 2 + p) x = 0 conservation of momentum (1-10) 

E t + (v(E + p)) x = 0 conservation of energy (1-11) 

The resulting system of three conservation laws gives the Euler equations of gas dynam- 
ics. Note that these equations involve another quantity, the pressure p, which must be 
specified as a given function of p, pv, and E in order that the fluxes are well defined func- 
tions of the conserved quantities alone. This additional equation is called the equation 
of state and depends on physical properties of the gas under study. 

If we introduce the vector u £ IR 3 as 



p(x,t) 


’ «1 " 

u(x,t) = 

p(x,t)v{x,t) 

= 

u 2 


E(x,t) 


. w 3 . 


( 1 . 12 ) 


then the system of equations (1.9), (1.10), (1.11) can be written simply as 

u t + f(u) x = 0 (1-13) 


where 



pv 


U 2 

f(u) = 

pv 2 -|- p 

= 

+ p(u) 


_ v(E + p) _ 


. U 2 (u 3 + p(u))/U! , 


(1.14) 


Again, the form (1.13) is the differential form of the conservation laws, which holds in 
the usual sense only where u is smooth. More generally, the integral form for a system of 
m equations says that 


/ u(x,t)dx = f(u(xi,t)) - f(u(x 2 ,t)) 
dt Jx i 


(1.15) 
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for all * 1 , x 2 , t. Equivalently, integrating from f x to t 2 gives 

u(x,t 2 )dx = j u(x,ti)dx ( 1 . 16 ) 

+ / f(u(x u t))dt- f 3 f(u(x 2 ,t))dt 
J h Jt x 

for all Xi, x 2 , and < 2 . These integral forms of the conservation law are the fundamental 
physical conservation laws. It is important to note that these equations make sense even 
if the function u(x,t) is discontinuous, while transforming to the differential form (1.13) 
is valid only when u is smooth. 

1*1.2 The shock tube problem 

A simple example that illustrates the interesting behavior of solutions to conservation laws 
is the “shock tube problem” of gas dynamics. The physical set-up is a tube filled with 
gas, initially divided by a membrane into two sections. The gas has a higher density and 
pressure in one half of the tube than in the other half, with zero velocity everywhere. At 
time t = 0, the membrane is suddenly removed or broken, and the gas allowed to flow. 
We expect a net motion in the direction of lower pressure. Assuming the flow is uniform 
across the tube, there is variation in only one direction and the one-dimensional Euler 
equations apply. 

The structure of this flow turns out to be very interesting, involving three distinct waves 
separating regions in which the state variables are constant. Across two of these waves 
there are discontinuities in some of the state variables. A shock wave propagates into the 
region of lower pressure, across which the density and pressure jump to higher values and 
all of the state variables are discontinuous. This is followed by a contact discontinuity, 
across which the density is again discontinuous but the velocity and pressure are constant. 
The third wave moves in the opposite direction and has a very different structure: all of 
the state variables are continuous and there is a smooth transition. This wave is called a 
rarefaction wave since the density of the gas decreases (the gas is rarefied) as this wave 
passes through. 

If we put the initial discontinuity at x = 0, then the resulting solution u{x y t) is 
a “similarity solution” in the variable x/t , meaning that u(x y t) can be expressed as a 
function of x/t alone, say t*(x,*) = w(x/t). It follows that u(x,t) = u(ax y at) for any 
Of > 0, so the solution at two different times t and at look the same if we rescale the 
x-axis. This also means that the waves move at constant speed and the solution ti(x, t) is 
constant along any ray x/t = constant in the x-t plane. 

Figure 1.1 shows a typical solution as a function of x/t . We can view this as a plot of 
the solution as a function of x at time t = 1, for example. The structure of the solution 
in the x-t plane is also shown. 
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density velocity 



pressure wave structure in x-t plane 



Figure 1.1. Solution to a shock tube problem for the one-dimensional Euler equations. 
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In a real experimental shock tube, the state variables would not be discontinuous 
across the shock wave or contact discontinuity because of effects such as viscosity and heat 
conduction. These are ignored in the Euler equations. If we include these effects, using the 
full Navier-Stokes equations, then the solution of the partial differential equations would 
also be smooth. However, these smooth solutions would be nearly discontinuous, in the 
sense that the rise in density would occur over a distance that is microscopic compared to 
the natural length scale of the shock tube. If we plotted the smooth solutions they would 
look indistinguishable from the discontinuous plots shown in Figure 1.1. For this reason 
we would like to ignore these viscous terms altogether and work with the simpler Euler 
equations. 

The shock tube problem is a special case of what is known mathematically as the 
Riemann problem. In general, the Riemann problem consists of the conservation law 
Ut + f(u)r = 0 together with the special initial data 



x < 0 
x > 0. 


(1.17) 


This data consists of two constant states separated by a discontinuity. It turns out that 
this problem is relatively easy to solve in general, and that the resulting wave structure 
reveals a lot about the structure of solutions more generally. Many numerical methods are 
based on solving Riemann problems, including the methods to be discussed in the next 
lecture. For this reason, on of the main goals of the present lecture is to understand the 
solution of this Riemann problem for general conservation laws. 


1.1.3 The equation of state 

In the Euler equations, the total energy E is often decomposed as 

E = ^pv 2 + pe. (1-18) 

The first term here is the kinetic energy, while pe is the internal energy. The variable e, 
internal energy per unit mass, is called the specific internal energy. Internal energy 
includes rotational and vibrational energy and possibly other forms of energy in more 
complicated situations. In the Euler equations we assume that the gas is in chemical and 
thermodynamic equilibrium and that the internal energy is a known function of pressure 
and density: 

« = e(p,/>). (1.19) 

This is the “equation of state” for the gas, which depends on the particular gas under 
study. 

For an ideal gas, internal energy is a function of temperature alone, e = c(T), and T 
is related to p and p by the ideal gas law, 


p = TZpT 


( 1 . 20 ) 
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where 1Z is a constant. To good approximation, the internal energy is simply proportional 
to the temperature, 

e = c v T , (1.21) 

where c v is a constant known as the specific heat at constant volume. Such gases are 
called polytropic. If energy is added to a fixed quantity of gas, and the volume is held 
constant, then the change in energy and change in temperature are related via 


de = c v dT. 


( 1 . 22 ) 


On the other hand, if the gas is allowed to expand while the energy is added, and 
pressure is held constant instead, not all of the energy goes into increasing the internal 
energy. The work done in expanding the volume 1/p by d{l/p) is pd(l/p) and we obtain 
another relation 

de + pd(l/p) = CpdT (1-23) 


or 

d(e + p/p) = c p dT 

where c p is the specific heat at constant pressure. The quantity 


h = e + p/p 


(1.24) 


(1.25) 


is called the enthalpy. For a polytropic gas, c p is also assumed to be constant so that 
(1.24) yields 

h = CpT. (1.26) 


Note that by the ideal gas law, 

c p - c v = U. (1.27) 

The equation of state for an polytropic gas turns out to depend only on the ratio of 
specific heats, usually denoted by 


7 — Cp/c v . 


(1.28) 


Note that T = p/TZp so that 


_ cT _ (Mp = __E 

~ v \n) P ( 7 - 


i )p 


(1.29) 


by (1.27) and (1.28). Using this in (1.18) gives the common form of the equation of state 
for a polytropic gas: 


E = 


+ -pv* 
7 - 1 ^ 2 ^ 


(1.30) 
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1.1.4 Isothermal flow 

Suppose we consider the flow of gas in a tube that is immersed in a bath at a constant 
temperature T, and assume that this bath maintains a constant temperature within the 
gas. Then the ideal gas law (1.20) reduces to 

p = a 2 p (1.31) 


where a 2 = TIT is a constant and a is the sound speed. Note that maintaining this constant 
temperature requires heat flux through the wall of the tube, and so energy is no longer 
conserved in the tube. But mass and momentum are still conserved and these equations, 
together with the equation of state (1.31), lead to the isothermal equations, 



pv 

pv 2 + a 2 p 


- 0 . 


X 


(1.32) 


This system of two equations is a particularly nice example to use in illustrating the 
nonlinear theory. The algebra is relatively simple and yet the behavior is analogous to 
what is seen for the full Euler equations. 


1.2 Scalar equations 

We begin our study of conservation laws by considering the scalar case. Many of the 
difficulties encountered with systems of equations are already encountered here, and a 
good understanding of the scalar equation is required before proceeding. 

1.2.1 The linear advection equation 

We first consider the linear advection equation, 


u t + au x = 0. (1.33) 

The Cauchy problem is defined by this equation on the domain — oo < x < oo, t > 0 
together with initial conditions 

u(x,0) = «o(a:)- (1-34) 

The solution to this problem is simply 

u(x,t) = uo(x - at) (1.35) 

for t > 0, as can be easily verified. As time evolves, the initial data simply propagates 
unchanged to the right (if a > 0) or left (if a < 0) with velocity a. The solution u(x,t) 
is constant along each ray x - at = x 0 , which are known as the characteristics of the 
equation. (See Fig. 1.2 for the case a > 0.) 
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Note that the characteristics are curves in the x-t plane satisfying the ordinary differ- 
ential equations x'(t ) = a, x(0) = x 0 . If we differentiate u (x,t) along one of these curves 
to find the rate of change of u along the characteristic, we find that 

= ut + au x (1.36) 

= 0 , 

confirming that u is constant along these characteristics. 

1.2.2 Domain of dependence 

Note that solutions to the linear advection equation (1.33) have the following property: 
the solution u(x, t ) at any point (x, t) depends on the initial data ti 0 only at a single point, 
namely the point xo such that (x, t) lies on the characteristic through io. We could change 
the initial data at any points other than x 0 without affecting the solution u(x,t). The set 
V(x,t) = {x 0 } is called the domain of dependence of the point (x,t). Here this domain 
consists of a single point. For a system of equations this domain is typically an interval, 
but a fundamental fact about hyperbolic equations is that it is always a bounded interval. 
The solution at ( x,t) is determined by the initial data within some finite distance of the 
point x. The size of this set usually increases with t, but we have a bound of the form 
V C {x : \x - x\ < a max t} for some value a max . Conversely, initial data at any given 
point x 0 can influence the solution only within some cone {x : |x - x 0 | < a m »xf} of the 
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0?,t) 




Figure l.S. Domain of dependence and range of influence. 


x-t plane. This region is called the range of influence of the point i 0 . See Figure 1.3 
for an illustration. We summarize this by saying that hyperbolic equations have finite 
propagation speed; information can travel with speed at most This has important 

consequences in developing numerical methods. 


1.2.3 Nonsmooth data 

In computing (1.36), we have assumed differentiability of u(x,t). However, from our 
observation that the solution along a characteristic curve depends only on the one value 
uo(zo)t it is clear that spatial smoothness is not required for this construction of u(x,t) 
from uo(z). We can thus define a “solution” to the PDE even if uq(x) is not a smooth 
function. Note that if u 0 (x) has a singularity at some point *o (a discontinuity in uq 
or some derivative), then the resulting u(x,t) will have a singularity of the same order 
along the characteristic curve through Xo, but will remain smooth along characteristics 
through smooth portions of the data. This is a fundamental property of linear hyperbolic 
equations: singularities propagate only along characteristics. 

If uq is nondifferentiable at some point then u(z, t ) is no longer a classical solution of 
the differential equation everywhere. However, this function does satisfy the integral form 
of the conservation law, which continues to make sense for nonsmooth u. Recall that the 
integral form is more fundamental physically than the differential equation, which was 
derived from the integral form under the additional assumption of smoothness. It thus 
makes perfect sense to accept this generalized solution as solving the conservation law. 

Other approaches can also be taken to defining this generalized solution, which extend 
better to the study of nonlinear equations where we can no longer simply integrate along 
characteristics. 

One possibility is to approximate the nonsmooth data uo(x) by a sequence of smooth 
functions Uo(x), with 


||«o - “olli < f 
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as e —► 0. Here || * ||i is the 1-norm, defined by 


IMIi = / K*)!***- 

*/— oo 


(1.37) 


For the linear equation we know that the PDE together with the smooth data tt£ has a 
smooth classical solution u ( (x,t) for all t > 0. We can now define the generalized solution 
by taking the limit of as c — ► 0. For example, the constant coefficient 

problem (1.33) has classical smooth solutions 


u { (x,t) = Uq(x - at) 

and clearly at each time t the 1-norm limit exists and satisfies 


u(x, t ) = lim Uq(x - at) = u 0 (x - at) 

£~>0 

as expected. 

Unfortunately, this approach of smoothing the initial data will not work for nonlinear 
problems. As we will see, solutions to the nonlinear problem can develop singularities even 
if the initial data is smooth, and so there is no guarantee that classical solutions with data 
Uq(x) will exist. 

A better approach, which does generalize to nonlinear equations, is to leave the initial 
data alone but modify the PDE by adding a small diffusive term. Recall that the Euler 
equations, for example, arise from the Navier-Stokes equations by ignoring the diffusive 
effects of viscosity and heat conduction. Analogously, the advection equation (1.33) can 
be considered as an approximation to the advection-diffusion equation 

u t + au x = eu xx (1.38) 

for c very small. If we now let u l (x, t) denote the solution of (1.38) with data u 0 (x), then 
u c is smooth even if Uo(x) is not smooth since (1.38) is a parabolic equation. We can 
again take the limit of u c (x,t) as e — ► 0, and will obtain the same generalized solution 
u(x, t) as before. Clearly this is the correct notion of a generalized solution to the inviscid 
equations, and is often called the “vanishing viscosity” solution. 


1.2.4 Burgers’ equation 

Now consider the nonlinear scalar equation 

«« + /(«)* = 0 (1.39) 

where f(u) is a nonlinear function of u. We will assume for the most part that /(u) is a 
convex function, f"(u) > 0 for all u (or, equally well, / is concave with f"{u) < 0 Vu). The 
convexity assumption corresponds to a “genuine nonlinearity” assumption for systems of 
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T / /V? T T t T T T T 


u.(y,o^ 

Figure 1 . 4 . Characteristics and solution for Burgers’ equation ( small t). 

equations that holds in many important cases, such as the Euler equations. The nonconvex 
case is also important in some applications (e.g. oil reservoir simulation) but is more 
complicated mathematically. 

By far the most famous model problem in this field is Burgers* equation, in which 
/( u) = Ju 2 , so (1.39) becomes 

u t + uu x — 0. (1-40) 

Actually this should be called the “inviscid Burgers’ equation”, since the equation origi- 
nally studied by Burgers also includes a viscous term: 

Uf \-uu x = eu xx . (1-41) 

This is about the simplest model that includes the nonlinear and viscous effects of fluid 
dynamics. 

Consider the inviscid equation (1.40) with smooth initial data. For small time, a 
solution can be constructed by following characteristics. Note that (1.40) looks like an 
advection equation, but with the advection velocity u equal to the value of the advected 
quantity. The characteristics satisfy 



x'(t ) = u(x(t),t) 

and along each characteristic u is constant, since 

^u(x(t),t) = jju(x(t),t) + £u(z(t),t)x'(t) 
= u t + uu x 
= 0 . 


(1.42) 


(1.43) 


Moreover, since u is constant on each characteristic, the slope x'(t) is constant by (1.42) 
and so the characteristics are straight lines, determined by the initial data (see Fig. 1.4). 
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U. (X, o') 


Figure 1.5. Shock formation in Burgers’ equation. 


If the initial data is smooth, then this can be used to determine the solution u(x, t) for 
small enough t that characteristics do not cross: For each (x,t) we can solve the equation 

s = £ + «({, 0)t (1.44) 


for £ and then 


u(x,t) = u(£,0). 


(1.45) 


1.2.5 Shock formation 


For larger t the equation (1.44) may not have a unique solution. This happens when the 
characteristics cross, as will eventually happen if u x (x,0) is negative at any point. At the 
time 7), where the characteristics first cross, the function u(x,t) has an infinite slope — 
the wave “breaks” and a shock forms. Beyond this point there is no classical solution of 
the PDE, and the weak solution we wish to determine becomes discontinuous. 

Figure 1.5 shows an extreme example where the initial data is piecewise linear and 
many characteristics come together at once. More generally an infinite slope in the solution 
may appear first at just one point x, corresponding via (1.44) to the point £ where the 
slope of the initial data is most negative. This gives the “breaking time” 


T b = 


-1 

min Uq(x)" 


(1.46) 


For times t > Tb some of the characteristics have crossed and so there are points x 
where there are three characteristics leading back to t = 0. One can view the “solution” 
u at such a time as a triple- valued function (see Fig. 1.6). 
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Figure 1.6. Triple-valued solution to Burgers’ equation at time t > 7j>. 

This sort of solution makes sense in some contexts, for example a breaking wave on 
a sloping beach can be modeled by the shallow water equations and, for a while at least, 
does behave as seen in Fig. 1.6, with fluid depth a triple- valued function. 

However, in most physical situations this does not make sense. For example, the 
density of a gas cannot possibly be triple valued at a point. What happens instead at 
time T{,? 

We can determine the correct physical behavior by adopting the vanishing viscosity 
approach. The equation (1.40) is a model of (1.41) valid only for small c and smooth 
u. When it breaks down, we must return to (1.41). If the initial data is smooth and c 
very small, then before the wave begins to break the eu xx term is negligible compared 
to the other terms and the solutions to both PDEs look nearly identical. Figure 1.4, 
for example, would be essentially unchanged if we solved (1.41) with small t rather than 
(1.40). However, as the wave begins to break, the second derivative term u xx grows much 
faster than u x , and at some point the cu xx term is comparable to the other terms and 
begins to play a role. This term keeps the solution smooth for all time, preventing the 
breakdown of solutions that occurs for the hyperbolic problem. 

For very small values of e, the discontinuous solution at Tj, shown in Figure 1.5 would 
be replaced by a smooth continuous function as in Figure 1.7. As e -» 0 this becomes 
sharper and approaches the discontinuous solution of Figure 1.5. 

For times t > Tb , such as was shown in Figure 1.6, the viscous solution for e > 0 would 
continue to be smooth and single valued, with a shape similar to that shown in Figure 
1.7. The behavior as c — ► 0 is indicated in Figure 1.8. 

It is this vanishing viscosity solution that we hope to capture by solving the inviscid 
equation. 
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Figure 1.7. Solution to the viscous Burgers’ equation at time Tb for the data shown in 
Figure 1.5. 



Figure 1.8. Solution to the viscous Burgers’ equation for two different values of e. 

1.2.0 Weak solutions 

A natural way to define a generalized solution of the inviscid equation that does not require 
differentiability is to go back to the integral form of the conservation law, and say that 
is a generalized solution if (1.8) is satisfied for all x\, X 2 , t%, t?. 

There is another approach that results in a different integral formulation that is often 
more convenient to work with. This is a mathematical technique that can be applied more 
generally to rewrite a differential equation in a form where less smoothness is required 
to define a “solution”. The basic idea is to take the PDE, multiply by a smooth “test 
function”, integrate one or more times over some domain, and then use integration by 
parts to move derivatives off the function u and onto the smooth test function. The result 
is an equation involving fewer derivatives on u, and hence requiring less smoothness. 

In our case we will use test functions 4> € Cq(1R x 1R). Here Cq is the space of func- 
tion that are continuously differentiable with “compact support”. The latter requirement 
means that <f>(x,t ) is identically zero outside of some bounded set, and so the support of 
the function lies in a compact set. 

If we multiply u t -f f x = 0 by 4>(x, <) and then integrate over space and time, we obtain 
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Now integrate by parts, yielding 

foo r+oo roo 

/ / [<j> t u + <f> x f(u)] dx dt = - / ^(x,0)u(x,0)dx. (1-48) 

JO J—oo J — oo 

Note that nearly all the boundary terms which normally arise through integration by 
parts drop out due to the requirement that <f> have compact support, and hence vanishes 
at infinity. The remaining boundary term brings in the initial conditions of the PDE, 
which must still play a role in our weak formulation. 

Definition 1.1. The function u(x, t) is called a weak solution of the conservation law 
if (1.4S) holds for all functions <f> G Cq(IR x 1R). 

The vanishing viscosity generalized solution we defined above is a weak solution in 
the sense of (1.48), and so this definition includes the solution we are looking for. Unfor- 
tunately, weak solutions are often not unique, and so an additional problem is often to 
identify which weak solution is the physically correct vanishing viscosity solution. Again, 
one would like to avoid working with the viscous equation directly, but it turns out that 
there are other conditions one can impose on weak solutions that are easier to check and 
will also pick out the correct solution. These are usually called entropy conditions by 
analogy with the gas dynamics case, where a discontinuity is physically realistic only if 
the entropy of the gas increases as it crosses the shock. 


1.2.7 The Riemann Problem 


The conservation law together with piecewise constant data having a single discontinuity is 
known as the Riemann problem. As an example, consider Burgers’ equation u t + uu x = 0 
with piecewise constant initial data 



x < 0 
x > 0. 


(1.49) 


The form of the solution depends on the relation between tq and u r . 


Case I. ui > u r . 

In this case there is a unique weak solution, 


u 



Ui 

u r 


X < st 
X > st 


where 

S = (ui + u r )/ 2 


(1.50) 

(1.51) 


is the shock speed, the speed at which the discontinuity travels. A general expression 
for the shock speed will be derived below. Note that characteristics in each of the regions 
where u is constant go into the shock (see Fig. 1.9) as time advances. 
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Figure 1.9. Shock wave . 


Case II. u/ < u r . 

In this case there are infinitely many weak solutions. One of these is again (1.50), 
(1.51) in which the discontinuity propagates with speed s. Note that characteristics now 
go out of the shock (Fig. 1.10) and that this solution is not stable to perturbations. If the 
data is smeared out slightly, or if a small amount of viscosity is added to the equation, 
the solution changes completely. 

Another weak solution is the rarefaction wave 

{ ui x < uit 

x/t urf < x < u r t (1.52) 

u r x > u r t 

This solution is stable to perturbations and is in fact the vanishing viscosity generalized 
solution (Fig. 1.11). 

1.2.8 Shock speed 

The propagating shock solution (1.50) is a weak solution to Burgers’ equation only if the 
speed of propagation is given by (1.51). The same discontinuity propagating at a different 
speed would not be a weak solution. 

The speed of propagation can be determined by conservation. If M is large compared 
to si then by (1.15), u(x , t ) dx must increase at the rate 

d fM 

dt J M u ^ x '^ dx ~ /(«/)-/(« r) (1.53) 

= ^(«/ + Ur)(ti/ - U r ) 

for Burgers’ equation. On the other hand, the solution (1.50) clearly has 
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so that 


d 

— / u(x y t) dx = s(ut — u r ). 
dt J 


(1.55) 


Comparing (1.53) and (1.55) gives (1.51). 

More generally, for arbitrary flux function f(u) this same argument gives the follow- 
ing relation between the shock speed s and the states tq and ti r , called the Rankine- 
Hugoniot jump condition: 


f(u t ) - f{u r ) = *(U/ - U r ). 

For scalar problems this gives simply 

s _ /(«<) - fM _ [f] 

Ul — u r [u] 

where [•] indicates the jump in some quantity across the discontinuity. Note that any jump 
is allowed, provided the speed is related via (1.57). 

For systems of equations, ui — u r and f(u r ) — /(iq) are both vectors while s is still 
a scalar. Now we cannot always solve for s to make (1.56) hold. Instead, only certain 
jumps ui — u r are allowed, namely those for which the vectors /(tq) — f(u r ) and iq — u r 
are linearly dependent. 

Example 1.1. For a linear system with f(u) = Au , (1.56) gives 


(1.56) 

(1.57) 


A(ui - U r ) ~ S(u t ~ U r ), 


(1.58) 


i.e.y ui — u r must be an eigenvector of the matrix A and s is the associated eigenvalue. 
For a linear system, these eigenvalues are the characteristic speeds of the system. Thus 
discontinuities can propagate only along characteristics, a fact that we have already seen 
for the scalar case. 


1.2.9 Entropy conditions 

As demonstrated above, there are situations in which the weak solution is not unique and 
an additional condition is required to pick out the physically relevant vanishing viscosity 
solution. The condition which defines this solution is that it should be the limiting solution 
of the viscous equation as e — ► 0, but this is not easy to work with. We want to find simpler 
conditions. 

For scalar equations there is an obvious condition suggested by Figures 1.9 and 1.10. 
A shock should have characteristics going into the shock, as time advances. A propa- 
gating discontinuity with characteristics coming out of it, as in Figure 1.10, is unstable 
to perturbations. Either smearing out the initial profile a little, or adding some viscosity 
to the system, will cause this to be replaced by a rarefaction fan of characteristics, as in 
Figure 1.11. This is the simplest version of the entropy condition: 
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Entropy condition: (For scalar convex conservation laws.) A discontinuity propa- 
gating with speed s given by (1.56) satisfies the entropy condition if 


f'(ui) > s > f\u r ). 


(1.59) 


Note that f\u) is the characteristic speed. For convex /, the Rankine-Hugoniot speed s 
from (1.57) must lie between /'(u/) and /'(u r ), so (1.59) reduces to simply the requirement 
that f\ui) > /'(tz r ), which again by convexity requires u\ > u r . 

For nonconvex fluxes, or systems of equation, more complicated entropy conditions are 
often used. We will not pursue these here. 

1.3 Linear Hyperbolic Systems 

We now begin to investigate systems of equations. We start with constant coefficient 
linear systems. Here we can solve the equations explicitly by transforming to characteristic 
variables. We will also obtain explicit solutions of the Riemann problem and introduce a 
“phase space” interpretation that will be very useful in our study of nonlinear systems. 
Consider the linear system 


u t + Au x = 0 (1.60) 

u(a;,0) = « 0 (*) 

where u : IR x IR — + IR”* and A £ IR mxm is a constant matrix. This is a system of 
conservation laws with the flux function f(u) = Au. This system is called hyperbolic if 
A is diagonalizable with real eigenvalues, so that we can decompose 

A = RAR~ 1 (1.61) 

where A = diag(Ai, A 2 , . . . , A m ) is a diagonal matrix of eigenvalues and R = [rj |r 2 | -|r m ] 

is the matrix of right eigenvectors. Note that AR = RA, i.e., 

At p - \ p r p for p = 1, 2, . . ., m. (1.62) 

The system is called strictly hyperbolic if the eigenvalues are distinct. We will always 
make this assumption as well. 

1.3.1 Characteristic variables 

We can solve (1.60) by changing to the “characteristic variables” 


v = R 1 u. 


(1.63) 
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Multiplying (1.60) by R 1 and using (1.61) gives 

R-'ut + AR-'us = 0 (1.64) 

or, since iZ -1 is constant, 

v t + Av x = 0. (1.65) 

Since A is diagonal, this decouples into m independent scalar equations 

(v p )t + Ap(u p ) x = 0, p = 1, 2, ..., m. (1.66) 

Each of these is a constant coefficient linear advection equation, with solution 

v p (x,t) = v p (x - Ap<,0). (1-67) 

Since v = R -1 u, the initial data for v p is simply the pth component of the vector 

u(a:,0) = R~ 1 uo(x). (1.68) 

The solution to the original system is finally recovered via (1.63): 

u(x,l) = Rv(x,t). (1.69) 

Note that the value v p (x,t) is the coefficient of r p in an eigenvector expansion of the 
vector u(a:,<), i.e., (1.69) can be written out as 

m 

u(x,t) = ^ v p (x, t)r p . (1-70) 

p= i 

Combining this with the solutions (1.67) of the decoupled scalar equations gives 

m 

«(x, t) = v p (x - A p £,0)r p . (1*71) 

p= i 

Note that depends only on the initial data at the m points x — \ p t, so the domain 

of dependence is given by V(x , t) = {x — x — p = 1, 2, . . . , m}. 

The curves x = x 0 + A p t satisfying x\t) = A p are the “characteristics of the pth family”, 
or simply “p-charactcristics”. These are straight lines in the case of a constant coefficient 
system. Note that for a strictly hyperbolic system, m distinct characteristic curves pass 
through each point in the x-t plane. The coefficient v p (x,t) of the eigenvector r p in the 
eigenvector expansion (1.70) of u(x,t) is constant along any p-characteristic. 
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1.3*2 The Riemann Problem 


For the constant coefficient linear system, the Riemann problem can be explicitly solved. 
We will see shortly that the solution to a nonlinear Riemann problem has a simple structure 
which is quite similar to the structure of this linear solution, and so it is worthwhile 
studying the linear case in some detail. 

The Riemann problem consists of the equation u t + Au x = 0 together with piecewise 
constant initial data of the form 


u(x, 0) 


{ ui x < 0 
u r x > 0 


(1.72) 


Recall that the general solution to the linear problem is given by (1.71). For the 
Riemann problem we can simplify the notation if we decompose ui and u r as 


ui 


= Y, a p 

p~\ 


P=1 


Then 


Up(x,°) 


{ a v x < 0 

i % x > o 


and so 

vp{x, t) 

If we let P(x,t) be the maximum value of p for which x — A p t > 0, then 


a p if a: — A p t < 0 
(3 P if x - A p t > 0. 


(1.73) 


(1.74) 

(1.75) 


P{x,t) m 

u(x,t)= J2 (3 p r p + ot p r p . 

V - 1 p=f > (x,*) + 1 


(1.76) 


The determination of u(x,t) at a given point is illustrated in Figure 1.12. In the case 
shown, V\ = /?i while V 2 = ot 2 and U 3 = 03 . The solution at the point illustrated is thus 


u{x,t) = (3\r\ + a 2 r 2 + a 3 r 3 . (1-77) 

Note that the solution is the same at any point in the wedge between the x f = Ai and 
x f = A 2 characteristics. As we cross the pth characteristic, the value of x — A p t passes 
through 0 and the corresponding v p jumps from a p to f3 p . The other coefficients u, ( i ^ j ) 
remain constant. 

The solution is constant in each of the wedges as shown in Figure 1.13. Across the pth 
characteristic the solution jumps with the jump given by 

[«] = (J3 p - a p )r p . (1.78) 

Note that these jumps satisfy the Rankine-Hugoniot conditions (1.56), since f(u ) = Au 
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Figure 1.13. Values of solution u in each wedge of x-t plane. 
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leads to 


[/] = A[u] 

= (ftp - a p )Ar v 
= A p [u] 

and A p is precisely the speed of propagation of this jump. The solution u(x,t) in (1.76) 
can alternatively be written in terms of these jumps as 

u(x,t) = u/+ (ftp — a p )r p 

X p <x/t 

= «r- Yi (^“ a p) r P 

A p <x/t 

It might happen that the initial jump u r — u\ is already an eigenvector of A, if u r — ui = 
(/?,■ — a t ) r » f° r some i . In this case a p = j3 p for p ^ i . Then this discontinuity simply 
propagates with speed A,, and the other characteristics carry jumps of zero strength. 

In general this is not the case, however, and the jump u r — u\ cannot propagate as 
a single discontinuity with any speed without violating the Rankine-Hugoniot condition. 
We can view “solving the Riemann problem” as finding a way to split up the jump u r — U{ 
into a sum of jumps 


(1.79) 

(1.80) 


U r Ul — (/?i Ou)f*l 4" * * * 4" (/^m (1.81) 

each of which can propagate at an appropriate speed A^ with the Rankine-Hugoniot con- 
dition satisfied. 

For nonlinear systems we solve the Riemann problem in much the same way: The 
jump u r — will usually not have the property that [/] is a scalar multiple of [u], but we 
can attempt to find a way to split this jump up as a sum of jumps, across each of which 
this property does hold. (Although life is complicated by the fact that we may need to 
introduce rarefaction waves as well as shocks.) In studying the solution of the Riemann 
problem, the jump in the pth family, traveling at constant speed A p , is often called the 
p-wave. 

1.3.3 The phase plane 

For systems of two equations, it is illuminating to view this splitting in the phase plane. 
This is simply the ui~u 2 plane, where u = (ui,u 2 ). Each vector u{x,t) is represented by 
a point in this plane. In particular, m and u T are points in this plane and a discontinuity 
with left and right states u/ and u T can propagate as a single discontinuity only if u r - u\ 
is an eigenvector of A, which means that the line segment from u\ to u r must be parallel 
to the eigenvector 7*1 or r 2 . Figure 1.14 shows an example. For the state u\ illustrated 
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Figure l.l\. The Hugoniot locus of the state u t consists of all states that differ from u t by 
a scalar multiple of or r 2 . 


there, the jump from ui to u r can propagate as a single discontinuity if and only if u r lies 
on one of the two lines drawn through ui in the direction rj and r 2 . These lines give the 
locus of all points that can be connected to ui by a 1-wave or a 2- wave. This set of states 
is called the Hugoniot locus. We will see that there is a direct generalization of this to 
nonlinear systems in the next chapter. 

Similarly, there is a Hugoniot locus through any point u T that gives the set of all points 
u\ that can be connected to u r by an elementary p-wave. These curves are again in the 
directions ri and r 2 . 

For a general Riemann problem with arbitrary u\ and u r , the solution consists of two 
discontinuities travelling with speeds Ai and A 2 , with a new constant state in between that 
we will call u m . By the discussion above, 

Mm = An + a 2 r 2 (1.82) 

so that u m — ui = (/?i — a\)ri and u T — u m — (/? 2 — a 2 )r 2 . The location of u m in the phase 
plane must be where the 1-wave locus through u\ intersects the 2- wave locus through u r . 
This is illustrated in Figure 1.15a,. 

Note that if we interchange u T and u\ in this picture, the location of u m changes as 
illustrated in Figure 1.15b. In each case we travel from u\ to u T by first going in the 
direction r\ and then in the direction r 2 . This is required by the fact that A] < A 2 since 
clearly the jump between u\ and u m must travel slower than the jump between u m and u r 
if we are to obtain a single- valued solution. 
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Figure LI 5. The new state u m arising in the solution to the Riemann problem for two 
different choices of ui and u r . 

For systems with more than two equations, the same interpretation is possible but 
becomes harder to draw since the phase space is now m dimensional. Since the m eigen- 
vectors r p are linearly independent, we can decompose any jump u r — into the sum of 
jumps in these directions, obtaining a piecewise linear path from u\ to u r in m-dimensional 
space. 


1.4 Nonlinear Systems 

Now consider a nonlinear system of conservation laws 


tit + /(«)* = 0, (1.83) 

where u : IR x IR — » IR m and / : IR m — ► IR m . This can be written in the quasilinear form 

u t + A{u)u x = 0 (1-84) 

where A(u) — f'(u ) is the m x m Jacobian matrix. Again the system is hyperbolic if 
A(u) is diagonalizable with real eigenvalues for all values of u, at least in some range where 
the solution is known to lie, and strictly hyperbolic if the eigenvalues are distinct for 
all u. We will assume that this system is strictly hyperbolic, and order the eigenvalues of 
A in increasing order, 


Ai < A 2 < • • * < A m 


(1.85) 
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Since the eigenvalues are distinct, the eigenvectors are linearly dependent. We choose 
a particular basis for these eigenvectors, {r p (u)}^_ 1 , usually chosen to be normalized in 
some manner, e.g. ||r p (u)|| = 1. 

In the previous section we constructed the solution to the general Riemann problem 
for a linear hyperbolic system of conservation laws. Our next goal is to perform a similar 
construction for the nonlinear Riemann problem. In the linear case the solution consists 
of m waves, which are simply discontinuities traveling at the characteristic velocities of 
the linear system. In the nonlinear case our experience with the scalar equation leads 
us to expect more possibilities. In particular, the physically relevant vanishing viscosity 
solution may contain rarefaction waves as well as discontinuities. We will first ignore the 
entropy condition and ask a simpler question: is it possible to construct a weak solution 
of the Riemann problem consisting only of m discontinuities propagating with constant 

speeds < s 2 < •«• < s m ? As we will see, the answer is yes for ||u/ - u r || sufficiently 
small. 

1.4.1 The Hugoniot locus 

Recall that if a discontinuity propagating with speed s has constant values u and u on 
either side of the discontinuity, then the Rankine-Hugoniot jump condition must hold, 

f{u)~ f(u) = s(u-u). (1.86) 

Now suppose we fix the point u (E lR m and attempt to determine the set of all points u 
which can be connected to u by a discontinuity satisfying (1.86) for some s. This gives a 
system of m equations in m + 1 unknowns: the m components of u, and s. This leads us 
to expect one parameter families of solutions. 

We know that in the linear case there are indeed m such families for any u. In the 
pth family the jump u - u is some scalar multiple of r p , the pth eigenvector of A. We can 
parameterize these families of solutions using this scalar multiple, say £, and we obtain 
the following solution curves: 


u p (ti;u) = « + £r p 

5 p(£i **) = 

for p = 1, 2, . . . , m. Note that u p (0; u) = u for each p and so through the point u in phase 
space there are m curves (straight lines in fact) of possible solutions. This is illustrated 
in Figure 1.15 for the case m = 2. The two lines through each point are the states that 
can be connected by a discontinuity with jump proportional to iq or r 2 . 

In the nonlinear case we also obtain m curves through any point u, one for each 
characteristic family. We again parameterize these curves by u p (£;it) with w p (0;?i) = u 
and let s p (f;u) denote the corresponding speed. To simplify notation, we will frequently 
write these as simply u p (£), s p (£) when the point u is clearly understood. 
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The Rankine-Hugoniot condition gives 

/(«p(0) - /(«) = 5 p(0(«p(0 - «)• (1-87) 

Differentiating this expression with respect to £ and setting £ = 0 gives 

/'(«)Up(0) = Sp(0)w'(0) (1.88) 

so that Wp(0) must be a scalar multiple of the eigenvector r p (u) of /'(u), while s p ( 0) = 
Ap(ti). The curve u p (£) is thus tangent to r p (u ) at the point ti. We have already observed 
this, by a slightly different argument, in discussing weak shocks in Chapter 1.3. For a 
system of m = 2 equations, this is easily illustrated in the phase plane. An example for 
the isothermal equations of gas dynamics is discussed below, see Figure 1.16. 

For smooth /, it can be shown using the implicit function theorem that these solution 
curves exist locally in a neighborhood of u, and that the functions u p and s p are smooth. 
See Lax[42] or Smoller[66] for details. These curves are called Hugoniot curves. The set of 
all points on these curves is often collectively called the Hugoniot locus for the point u. 
If u p lies on the pth Hugoniot curve through u, then we say that u and u p are connected 
by a p-shock. 

Example 1.2. The isothermal equations of gas dynamics (1.32) provide a relatively 
simple example of the nonlinear theory. 

If we let m represent the momentum, m = pv , then the system becomes 


pt + m x - 0 
m t + (m 2 /p + a 2 p) x = 0 


or u t + }{u) x = 0 where u = (p, m). 
The Jacobian matrix is 


0 1 
a 2 - m 2 j p 2 2 m/p 


The eigenvalues are 
with eigenvectors 


/'(«) : 

\i(it) = m/p-a, A 2 (u) = m/p + a 


n(u) = 


l 

mj p — a 


r 2 (u) = 


1 

m/p + a 


(1.89) 


(1.90) 

(1.91) 

(1.92) 


Now let’s fix a state u and determine the set of states u that can be connected by a 
discontinuity. The Rankine-Hugoniot condition (1.86) becomes, for this system, 


m — m = s(p — p) 

(m 2 / p + a 2 p) - (m 2 / p + a 2 p) = s(m — m). 


(1.93) 
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This gives two equations in the three unknowns p, m, and s. These equations can be 
easily solved for rh and s in terms of p, giving 


m — pm/p ± a 



(1.94) 


s = rh/p ± a\J p/p. (1.95) 

The ± signs in these equations give two solutions, one for each family. Since m and s can 
be expressed in terms of p, we can easily parameterize these curves by taking, for example, 


Pp( 6«) = P( 1 + 0. p=l, 2. 


(1.96) 


We then have 


«i(6«) = « + £ 


P 

m - ap\J 1 + \ 


and 


u 2 (£;u) = u + ( 


P 

m + ap\J 1 + % ’ 


«i(C; «) = ™fp - “n/ 1 + 1- (1-97) 


«2(^; *) = m/p + a \fi + I- (1.98) 


The choice of signs for each family is determined by the behavior as £ — » 0, where the 
relation (1.88) must hold. It is easy to check that with the above choice we have 


d 

di 


ilp(0;u) 


s p (0;u) 


pr p {u) a r p (u), 
A p (m), 


as expected. 

Notice that real- valued solutions exist only for £ > —p and that u p {—p\u) = (0,0) for 
p = 1,2 and any u. Thus, each Hugoniot locus terminates at the origin (the vacuum state , 
since p = 0). There are no states with p < 0 that can be connected to u by a propagating 
discontinuity. The curves u p (f) are illustrated in Figure 1.16a for one particular choice 
of u and a — 1. Figure 1.16b shows how these curves vary with u. In this case we see 
« p (f;6) for u - (p, 0), p = 1, 2, 3, 4. 


1*4.2 Solution of the Riemann problem 

Now suppose that we wish to solve the Riemann problem with left and right states ui and 
u r (and recall that we are ignoring the entropy condition at this point). Just as in the 
linear case, we can accomplish this by finding an intermediate state u m such that u\ and 
u m are connected by a discontinuity satisfying the Rankine-Hugoniot condition, and so 
are u m and u r . Graphically we accomplish this by drawing the Hugoniot locus for each 
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A 

m 



Figure 1.16 . a) Hugoniot locus for the state u = (1,1) mi the isothermal gas dynamics 
equations with a = 1, b) Variation of these curves for u = (p, 0), /5 = 1, 2, 3, 4. 



Figure 7.7 7. Construction of a weak solution to the Riemann problem with left and right 
states ui and u r . 
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of the states u\ and u r and looking for intersections. See Figure 1.17 for an example with 
the isothermal equations. 

In this example there are two points of intersection, labelled u m and but only u m 
gives a single- valued solution to the Riemann problem since we need the jump from «/ to 
u m to travel more slowly than the jump from u m to u r . This requires that u m be connected 
to ui by a 1-shock while u r is connected to u m by a 2-shock, due to our convention that 
Ai(u) < A 2 (u). 

The state u m can be found algebraically by using our explicit expressions for the 
Hugoniot locus. We wish to find a state (p m ,rn m ) which is connected to tq by a 1-shock 
and to u r by a 2-shock. Consequently equation (1.94) with the minus sign should hold for 
u — u m , it = tq, and the same equation with the plus sign should hold for u = u my u = u r . 
Equating the two resulting expressions for m m gives 

PmTfti/ pi ~ ayjp m / pi (Pm ~ Pi) ~ Pm m r / pr 4“ Pm/ Pr (Pm ~ Pr )• (1.99) 


Setting z = we obtain a quadratic equation for z y 


a 

VFr 


+ 



^)z-a( v ^:+ v ^7) = o. 
Pi / 


( 1 . 100 ) 


This equation ha^ a unique positive solution z , and p m = z 2 . We can then determine m m 
by evaluating either side of (1.99). 

More generally, for a system of m equations we can attempt to solve the Riemann 
problem by finding a sequence of states u\, «2> •••? 1 such that u\ is connected to 

u\ by a 1-shock, U\ is connected to U2 by a 2-shock, and so on, with u m _i connected 
to u r by an m-shock. If u\ and u r are sufficiently close together then this can always 
be achieved. Lax proved a stronger version of this in his fundamental paper [41]. (By 
considering rarefaction waves also, the entropy satisfying solution can be constructed in a 
similar manner.) 


1.4.3 Genuine nonlinearity 

In defining the Hugoniot locus above, we ignored the question of whether a given dis- 
continuity is physically relevant. The state u is in the Hugoniot locus of it if the jump 
satisfies the Rankine-Hugoniot condition, regardless of whether this jump could exist in a 
vanishing viscosity solution. We would now like to define an entropy condition that can be 
applied directly to a discontinuous weak solution to determine whether the jumps should 
be allowed. Lax[41] proposed a simple generalization of the entropy condition (1.59) to 
systems of equations that are genuinely nonlinear (a natural generalization of the convex 
scalar equation). The pth characteristic field is said to be genuinely nonlinear if 


VAp(u) • r p (tt) / 0 for all u, 


( 1 . 101 ) 
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where VA p (u) = ( dX p /du\ y . .., d\ p /du m ) is the gradient of A p (u). Note that in the 
scalar case, m = 1 and X\(u) = /'(u) while ri(u) = 1 for all u. The condition (1.101) 
reduces to the convexity requirement /"(«) ^ 0 Vu. This implies that the characteristic 
speed /'(u) is monotonically increasing or decreasing as u varies, and leads to a relatively 
simple solution of the Riemann problem. 

For a system of equations, (1.101) implies that A p (u) is monotonically increasing or 
decreasing as u varies along an integral curve of the vector field r p (u). An integral curve 
of r p (u) is a curve with the property that it is everywhere tangent to r p (u). This will 
be discussed in detail shortly, where we will see that through a rarefaction wave u varies 
along an integral curve. Since monotonicity of the propagation speed A p is clearly required 
through a rarefaction wave, genuine nonlinearity is a natural assumption. 


1.4.4 The Lax entropy condition 

For a genuinely nonlinear field, Lax’s entropy condition says that a jump in the pth field 
(from tij to u r , say) is admissible only if 

Ap(tij) > 5 > A p (u r ) (1.102) 


where s is again the shock speed. Characteristics in the pth family disappear into the 
shock as time advances, just as in the scalar case. 

EXAMPLE 1.3. For isothermal gas dynamics we can easily verify that both fields are 
genuinely nonlinear. Since A p = m/p ± a, we compute that in each case 


VA p (u) 


Using (1.92), we compute that 


-m/p 2 

Up 


P= 1* 2. 


VAj(tt) ■ rj(u) = -a/p 
VA 2 (u) • r 2 (u) = a/p. 


(1.103) 


(1.104) 

(1.105) 


These quantities are both nonzero for all u . 

Now suppose ui and u r are connected by a 1-shock. Then u/ lies in the Hugoniot locus 
of u r and also u r lies in the Hugoniot locus of We can thus evaluate the shock speed 
s using (1.95) (with the minus sign, since the jump is a 1-shock) in two different ways, 
obtaining 

s=—~ aJpr/pi = — - aJ pi/ p T . (1.106) 

Pi V Pr V 

Since Ai(u) = m/p - a, the entropy condition (1.102) becomes 

mi j — m T / — m T 

a > aJp r /pi = aJpr/pi > a 

Pi PI V Pr * Pr 


(1.107) 
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and is clearly satisfied if and only if p r > /?/. 

Notice that since the fluid velocity is v = m/p, 1-shocks always travel slower than the 
fluid on either side, and so a given fluid particle passes through the shock from left to 
right (i.e. its state jumps from u\ to u T ). A consequence of the entropy condition is that 
the density of the gas must increase as it goes through the shock. This is also true more 
generally in the full Euler equations. The gas can only be compressed as the shock passes, 
not rarefied (rarefaction occurs, naturally enough, through a rarefaction wave rather than 
a shock). 

For 2-shocks the entropy condition requires 

— + a > — + ayffrjpi = — + a\J ~pjjpi > — + a (1.108) 

Pi PI v Pr v Pr 

which is now satisfied only if p T < pi. But note that 2-shocks travel faster than the fluid 
on either side, so that particles pass through the shock from right to left. So the entropy 
condition has the same physical interpretation as before: the density must jump from p r 
to a higher value pi as the gas goes through the shock. 

We can now reconsider the Hugoniot locus of a point u and retain only the points u 
that can be connected to u by an entropy-satisfying shock, discarding the entropy-violating 
shocks. In order to do this, we must first decide whether u is to lie to the left of the 
discontinuity or to the right. The entropy condition (1.102), unlike the Rankine-Hugoniot 
condition, is not symmetric in the two states. 

Figure 1.18a shows the set of states that can be connected to the right of a given 
state u by an entropy-satisfying shock. Figure 1.18b shows the set of states that can be 
connected to the left of the same state it . Note that the union of these curves gives the full 
Hugoniot locus. Each branch of the Hugoniot locus splits into two parts at ti; states on 
one side can only be connected to the left, states on the other side can only be connected 
to the right. 

1.4.5 Linear degeneracy 

The assumption of genuine nonlinearity is obviously violated for a constant coefficient 
linear system, in which A p (u) is constant and hence VA P = 0. More generally, for a 
nonlinear system it might happen that in one of the characteristic fields the eigenvalue 
A p (u) is constant along integral curves of this field, and hence 

VA p (u) • r p (u) = 0 Vu. (1-109) 

(Of course the value of A p (u) might vary from one integral curve to the next.) In this case 
we say that the pth field is linearly degenerate. This may seem rather unlikely, and not 
worth endowing with a special name, but in fact the Euler equations have this property. 
For this system of three equations, two of the fields are genuinely nonlinear while the third 
is linearly degenerate. 
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Figure 1.18. a) States u r that can be connected to u\ — u by an entropy-satisfying shock, 
b) States u/ that can be connected to u T = u by an entropy-satisfying shock . In each case 
the entropy-violating branches of the Hugoniot locus are shown as dashed lines. 

A discontinuity in a linearly degenerate field is called a contact discontinuity. This 
name again comes from gas dynamics. In a shock tube problem the gas initially on one 
side of the diaphragm never mixes with gas from the other side (in the inviscid Euler 
equations). As time evolves these two gases remain in contact along a ray in the x-t plane 
along which, in general, there is a jump in density. This is the contact discontinuity. 

1.4.0 Rarefaction Waves and Integral Curves 

We now turn our attention to rarefaction waves. All of the Riemann solutions considered 
so far have the following property: the solution is constant along all rays of the form 
x — £t. Consequently, the solution is a function of x/t alone, and is said to be a “similarity 
solution” of the PDE. A rarefaction wave solution to a system of equations also has this 
property and takes the form 

U{ x < t 

u(x,t)=< w(x/t) i\t < x < (1.110) 

k U T X > ( 2 1 

where w is a smooth function with w(£i) = uj and ^(£ 2 ) = u r . 

When does a system of equations have a solution of this form? As in the case of shocks, 
for arbitrary states ui and u r there might not be a solution of this form. But in general, 
starting at each point ui there are m curves consisting of points u r which can be connected 
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to Uj by a rarefaction wave. These turn out to be subsets of the integral curves of the 
vector fields r p (u). 

An integral curve for r p (u ) is a curve which has the property that the tangent to the 
curve at any point u lies in the direction r p (u). The existence of smooth curves of this 
form follows from smoothness of / and strict hyperbolicity, since r p (u) is then a smooth 
function of u. If u p (£) is a parameterization (for £ e IR) of an integral curve in the pth 
family, then the tangent vector is proportional to r p (u p (£ )) at each point, i.e. 

u p(0 = «(0 r p(« P (0) (1.111) 

where a(£) is some scalar factor. 

To see that rarefaction curves lie along integral curves, and to explicitly determine the 
function w(x/t) in (1.110), we differentiate u(x,t) = w{x/i) to obtain 

u t (x,t) = ~^w'{x/t) (1.112) 

u r (x,t) = ^w\x/t) (1.113) 

so that u t + f'{u)u x - 0 yields 

- -pw'(x/t) + jf'(w(x/t))w'(x/t) = 0. (1.114) 

Multiplying by t and rearranging gives 

= MO. (i.ii5) 

where £ = xjt. one possible solution of (1.115) is u/(£) = 0, i.e., w constant. Any constant 
function is a similarity solution of the conservation law, and indeed the rarefaction wave 
(1.110) takes this form for £ < ft and £ > £ 2 . In between, w is presumably smoothly 
varying and w / 0. Then (1.115) says that u/(£) must be proportional to some eigenvector 

Mt 0 (O)of/'(to(O), 

w> (0 = <x(t)r p (w(Q) (1.116) 

and hence the values w($) all lie along some integral curve of r p . In particular, the states 
u\ = w(£i) and u r — w(£ 2 ) both lie on the same integral curve. This is a necessary 
condition for the existence of a rarefaction wave connecting m/ and u T , but note that it is 
not sufficient. We need £ = x/t to be monotonically increasing as u;(^) moves from u / to 
u r along the integral curve; otherwise the function (1.110) is not single- valued. Note that 
our parameterization of the integral curve by £ is not at all arbitrary at this point, since 
(1.115) requires that £ be an eigenvalue of f(w( £)), 


e = MMO). 


( 1 . 117 ) 
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This particular parameterization is forced by our definition £ = x/t. Note that (LI 17) 
implies that w is constant along the ray x = A p (tn)tf, and hence each constant value of w 
propagates with speed A p (tn), just as in the scalar problem. 

By (1.117), monotonicity of £ is equivalent to monotonicity of A p (ti;) as w moves from 
ui to u T . From a given state u/ we can move along the integral curve only in the direction 
in which A p is increasing. If A p has a local maximum at «/ in the direction r p , then there 
are no rarefaction waves with left state «/. In the generic nonlinear case, there is a one 
parameter family of states that can be connected to u\ by a p-rarefaction - all those states 
lying on the integral curve of r p in the direction of increasing A p up to the next local 
maximum of A p . 

If the pth field is genuinely nonlinear then A p is monotonically varying along the entire 
integral curve. We need not worry about local maxima and we see that u/ and u r can 
always be connected by a rarefaction wave provided they lie on the same integral curve 
and 

A p (tq) < A p (ti r ). (1.118) 

If the pth field is linearly degenerate, then A p is constant on each integral curve and 
there are no possible rarefaction waves in this family. 

In order to explicitly determine the function w(f), we first determine the scale factor 
a(£) in (1.116) by differentiating (1.117) with respect to £. This gives 

1 = VA p (w(0) • u>'(0 

= «(0 VA pM0)- r P ( w (0) 


using (1.116), and hence 

“ (0 = VA p (,»(f)).r p («,(0)- 

Using this in (1.116) gives a system of ordinary differential equations for w(£): 


(1.119) 


w '(0 


r p( w (0) 

VA P MO)-rXO)’ 


6 < £ < 6 


( 1 . 120 ) 


with initial data 

w(6) = “I 

where = A p (u;) and £2 = \>(M r ). Note that the denominator in (1.120) is finite for 
f < £2 only if ^p is monotone between £1 and £ 2 - 

Example 1.4. We will construct 1-rarefactions for the isothermal equations. Using 
(1.92) and (1.104), the system of ODEs (1.120) takes the form 


p'(0 = ~p(0/ a i P(6) = Pi 

™'(0 = p(0 - m(£,) = m, 


( 1 . 121 ) 
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where = Ai(tq) = mi/pi - a. The first ODE is decoupled from the second and has 
solution 

p(0 = (1.122) 

The second ODE is then 

m'(£) = — m(£)/a, m(fj) = m/ (1.123) 


with solution 


m( 0 = (M^-6) + "*/)e" K " {,)/ ° 

= p,(t + a)e-«-M a . (1.124) 

From the solutions (p(£),m(£)) lt * s a ^ so use f u l to eliminate f and solve for m as a 
function of p. This gives explicit expressions for the integral curves in the phase plane. If 
we solve for £ in (1.122) and use this in (1.124) we obtain 

m{p) - pmi/pi - ap\og(p/pi). (1.125) 

We can construct 2- rarefactions in exactly the same manner, obtaining 

p(t) = pie«~^'°, (1.126) 

m(0 = pi(£ — a)e^-^ 1 ^ a . (1.127) 

and consequently 

m(p) = pm-i/pi + ap\og(p/pi). (1.128) 

For a given state u = u/ we can plot the set of all states u T which can be connected to 
uj by a rarefaction wave in either the first or second family. This is shown in Figure 1.19a 
for a particular choice of tq. Note that if we consider this same state u to be u T and now 
plot the set of all states u t that can be connected to u = u r by a rarefaction, we obtain 
a different picture as in Figure 1.19b. We must now have f decreasing as we move away 
from u T and so it is the opposite side of each integral curve that is now relevant. 

Note that these integral curves are very similar to the Ilugoniot locus, e.g\, Figure 1.18. 
Locally, near the point u, they must in fact be very similar. We know already that in the 
pth family each of these curves is tangent to r p (u) at u. Moreover, it can be shown that 
the curvature of both curves is the same (See Lax[42]). 

1.4.7 General solution of the Riemann problem 

We can combine Figures 1.18 and 1.19 to obtain a plot showing all states that can be 
connected to a given u by entropy-satisfying waves, either shocks or rarefactions. Again, 
the nature of this plot depends on whether u is to be the left state or right state, so we 



38 


1 Hyperbolic Conservation Laws 


2 -r&reW: 


t*) |(^ 

Figure 1.1 9. a) Set of states that can be connected to u\ = u by a rarefaction wave . b) 
Set of states that can be connected to u r = u by a rarefaction wave . In each case the full 
integral curves are shown as dashed lines . 

obtain two plots as shown in Figure 1.20. Here S p is used to denote the locus of states 
that can be connected by a p-shock and 7 Z p is the locus of states that can be connected 
by a p-rarefaction. Notice that the shock and rarefaction curves match up smoothly (with 
the same slope and curvature) at the point u. 

To solve the general Riemann problem between two different states u\ and u r , we 
simply superimpose the appropriate plots and look for the intersection u m of a 1-wave 
curve from u/ and a 2- wave curve from u r . An example for the isothermal equations is 
shown in Figure 1.21. This is the same example considered in Figure 1.17. We now see 
that the entropy-satisfying weak solution consists of a 1-shock from u\ to u m followed by 
a 2- rarefaction from u m to u T . 

To analytically determine the state u m , we must first determine whether each wave is 
a shock or rarefaction, and then use the appropriate expressions relating m and p along 
each curve to solve for the intersection. We have already seen how to do this for the case 
of two shocks, by solving the equation (1.99). If the solution consists of two rarefactions 
then the intermediate state must satisfy 

m m — PmlTll/ Pi l°S(/^m/P/) (1.129) 

rn m = p m m r j p r + aprn log(p m /p r ). (1.130) 

Equating the two right hand sides gives a single equation for p m alone, with solution 

(l (mi m r \\ 




(1.131) 
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Figure 1.20 . a) Set of states that can be connected to «/ by an entropy- satisfying 1-wave 
or 2-wave . b) Set of states that can be connected to u r . In each case f 1 Z p denotes 
p-rarefactions and S p denotes p-shocks. 



Figure 1.21. Construction of the entropy-satisfying weak solution to the Riemann problem 
with left and right states ui and u r . 
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We then obtain m m from either (1.129) or (1.130). 

If the solution consists of one shock and one rarefaction wave, as in Figure 1.21, then 
we must solve for p m and m m from the equations 

m„ = ^ - a .[^ - p,) (1.132) 

pi V Pi 

m m = + ap m log (p m l pi ) , 

pr 

for example, in the case of a 1-shock followed by a 2-rarefaction. In this case it is not 
possible to obtain a closed form solution ( p my m m ). Instead, it is necessary to solve these 
two equations by an iterative method such as Newton’s method. 

1*5 The Riemann problem for the Euler equations 

If we compute the Jacobian matrix f'(u) from (1.14), with the polytropic equation of state 
(1.30), we obtain 


/'(«) 


0 

-|(7 + l)i> 2 

-v(E + p)/p+ i (7 - l)t> 3 


1 

(3 — i)v 

(E + P)/P~( l~l)v 2 


0 ' 

(7-1) • 

1 V 


The eigenvalues are 


(1.133) 


A x(u) = v - c, A 2 (u) = v , A 3 (u) - v + c 

where c is the sound speed, 



(1.134) 

(1.135) 


1.5.1 Contact discontinuities 

Of particular note in these equations is the fact that the second characteristic field is 
linearly degenerate. It is easy to check from (1.133) that 


r 2 {u) 


1 


V 


L 2 V 


(1.136) 


is an eigenvector of /'(it) with eigenvalue A 2 (u) = v = (pv)/ p. Since 

' -v/p 

va 2 (w)= \jp 

o 


(1.137) 
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Figure 1.22. Typical solution to the Riemann problem for the Euler equations. 


we find that VA 2 ■ r 2 — 0. 

Since the second field is linearly degenerate, we can have neither rarefaction waves nor 
shocks in this field. Instead we have contact discontinuities, which are linear discontinuities 
that propagate with speed equal to the characteristic speed A 2 , which is simply the fluid 
velocity v. Across a contact discontinuity there is a jump in the density of the gas but 
the pressure and velocity are continuous. 

1.5.2 Solution to the Riemann problem 

The first and third characteristic fields are genuinely nonlinear and have behavior similar 
to the two characteristic fields in the isothermal equations. The solution to a Riemann 
problem typically has a contact discontinuity and two nonlinear waves, each of which 
might be either a shock or a rarefaction wave depending on u\ and u r . A typical solution 
is shown in Figure 1.22. 

Because v and p are constant across the contact discontinuity, it is often easier to work 
in the variables (p, u,p) rather than (p,pu, f?), although of course the jump conditions 
must be determined using the conserved variables. The resulting Hugoniot locus and 
integral curves can be transformed into (p, u,p) space. 

If the Riemann data is (p/, tv,p/) and (p r ,tv,p r ), then the two new constant states that 
appear in the Riemann solution will be denoted by u * = (p*,t>*,p*) and u * = (p*,v*,p*). 
(See Figure 1.22.) Note that across the 2- wave we know there is a jump only in density. 

Solution of the Riemann problem proceeds in principle just as before. Given the states 
ui and u r in the phase space, we need to determine the two intermediate states in such a 
way that u/ and u * are connected by a 1-wave, u* and u* are connected by a 2- wave, and 
finally u* and u r are connected by a 3-wave. 

This seems difficult, but we can take advantage of the fact that we know the 2-wave will 
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Figure 1.23 . Projection of shock and rarefaction curves onto the two-dimensional v-p 
plane y and determination of u* . 

be a contact discontinuity across which v and p are constant to make the problem much 
simpler. Instead of considering the full three dimensional (p, v,p) phase space, consider 
the v-p plane and project the integral curves and Hugoniot loci for the 1-waves and 3- 
waves onto this plane. In particular, project the locus of all states that can be connected 
to itf by a 1-wave (entropy satisfying shocks or rarefactions) onto this plane and also the 
locus of all states that can be connected to u r by a 3- wave. This gives Figure 1.23. 

We see in this example that we can go from ut (or actually, the projection of ui) to u* 
by a 1-rarefaction and from u* to u r by a 3-shock. The problem with this construction, of 
course, is that these curves are really curves in 3-space, and just because their projections 
intersect does not mean the original curves intersect. However, the curve 1Z\(ui) must go 
through some state u* = (p*, v*,p*) for some p * (so that it’s projection onto the v-p plane 
is (t>*,p*)). Similarly, the curve <? 3 (u r ) must pass through some state u* = (p*, v*,p*). 
But these two states differ only in p, and hence can be connected by a 2-wave (contact 
discontinuity). We have thus achieved our objective. Note that this technique depends on 
the fact that any jump in p is allowed across the contact discontinuity. 

In practice the calculation of u * can be reduced to a single nonlinear equation for 
p*, which is solved by an iterative method. Once p * is known, u* y pf and p* are easily 
determined. Godunov first proposed a numerical method based on the solution of Riemann 
problems and presented one such iterative method in his paper[29] (also described in §12.15 
of [59]). Chorin[16] describes an improvement of this method. More details on the solution 
of the Riemann problem can also be found in §81 of [22]. 



